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ABSTRACT 

Recent surveys of star forming regions have shown that most stars, and proba- 
bly all massive stars, are born in dense stellar clusters. The mechanism by which a 
molecular cloud fragments to form several hundred to thousands of individual stars 
has remained elusive. Here, we use a numerical simulation to follow the fragmentation 
of a turbulent molecular cloud and the subsequent formation and early evolution of a 
stellar cluster containing more than 400 stars. We show that the stellar cluster forms 
through the hierarchical fragmentation of a turbulent molecular cloud. This leads to 
the formation of many small subclusters which interact and merge to form the final 
stellar cluster. The hierarchical nature of the cluster formation has serious implica- 
tions in terms of the properties of the new-born stars. The higher number-density of 
stars in subclusters, compared to a more uniform distribution arising from a mono- 
lithic formation, results in closer and more frequent dynamical interactions. Such close 
interactions can truncate circumstellar discs, harden existing binaries, and potentially 
liberate a population of planets. We estimate that at least one-third of all stars, and 
most massive stars, suffer such disruptive interactions. 

Key words: stars: formation - stars: luminosity function, mass function - globular 
clusters and associations: general. 



1 INTRODUCTION 

The advent of large, efficient infrared detectors has resulted 
in a fundamental shift in our understanding of star for- 
mation. From the older viewpoint that star formation was 
something done independently and in isolation (Shu, Adams 
& Lizano 1987) , we now know that star formation is a group 
activity whereby some tens to thousands of stars form within 
a fraction of a parsec of each other (Clarke, Bonnell & Hil- 
lenbrand 2000). In such environments, the forming stars in- 
teract with each other on timescales comparable to their for- 
mation, resulting in a highly dynamical picture of star for- 
mation (Bate, Bonnell & Bromm 2003; Bonnell et al. 1997; 
Chapman et al. 1992). Infrared surveys of star forming re- 
gions have repeatedly shown that most stars form in clus- 
ters (Lada et al. 1991; Lada 1999; Clarke et al. 2000). This 
is found from both unbiased and pointed observations of 
molecular clouds and it is estimated that between 50 and 
95 per cent of stars form in clusters. Massive stars are even 
more likely to be found in young stellar clusters (Clarke 
et al. 2000). Testi et al.(1997) used pointed observations of 
Herbig AeBe stars and found a clear relation between the 
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star's mass and a surrounding stellar cluster. This implies a 
potential causal relationship between the cluster properties 
and the mass of the most massive star. 

The formation of stellar clusters has been an unsolved 
problem in astronomy due to the intrinsic difficulty of frag- 
menting a molecular cloud into a large number of stars, cou- 
pled with the numerical difficulties in following the subse- 
quent dynamical evolution. The spherical, nearly homoge- 
neous nature of young stellar clusters has generally been 
thought to imply that the preceding molecular cloud was 
itself smooth and nearly spherical (Goodwin 1998). Frag- 
menting such an object is extremely difficult as it requires 
the existence of self-gravitating clumps that have gas den- 
sities significantly higher than the mean gas density of the 
cloud (Bonnell 1999). 

In previous simulations involving smaller, lower-mass 
clouds, fragmentation into many bodies occurs most readily 
when the cloud contains filamentary structures which easily 
satisfy the above condition (Bastien et al. 1991; Inutsuka & 
Miyama 1997: Klessen, Burkert & Bate 1998; Bonnell 1999; 
Klessen 2000 & Burkert 2000). Filamentary structures oc- 
cupy a small fraction of the total volume of the cloud. Thus, 
the free-fall time of any self-gravitating perturbation in the 
structure is shorter than the overall dynamical time of the 



2 LA. Bonnell et al. 



system. Thus they can collapse to form fragments before col- 
liding together, forming binary and multiple systems (Bon- 
nell et al. 1991). A recent calculation (Bate et al. 2003) mod- 
elled the formation of « 50 low-mass stars and brown dwarfs 
from the fragmentation of a turbulent molecular cloud. In 
this simulation (see also Klessen et al. 1998; Klessen & Burk- 
ert 2000), the turbulence leads to filamentary structures 
which then fragment as stated above. Cloud-cloud collisions 
have also been shown to lead to sheet-like and then filamen- 
tary structures which subsequently fragment into multiple 
systems (Chapman et al. 1992). 

In this paper, we present the results from the first nu- 
merical simulation to follow the formation and evolution of a 
cluster containing more than 400 stars. The primary new in- 
sight from this simulation is that the cluster forms through a 
hierarchical process of many small sub-clusters which grow 
and merge through the subsequent dynamics to form the 
much larger final cluster. In Section 2 we detail the calcu- 
lation. In Section 3 we describe the evolution of the form- 
ing cluster while in section 4 we discuss the resultant mass 
distribution. Section 5 describes the implications of a hier- 
archical cluster formation process for stellar properties. Our 
conclusions are given in section 6. 



2 CALCULATIONS 

We use a high resolution Smoothed Particle Hydrodynamics 
(SPH) (Monaghan 1992) simulation to follow the fragmenta- 
tion of an initially uniform density molecular cloud contain- 
ing 1000A/q in a diameter of 1 pc and a gas temperature of 
10 K. SPH is a Lagrangian particle-based method that calcu- 
lates fluid properties by interpolation. Calculations of grav- 
itational forces are facilitated using a tree-structure (Benz 
et al. 1990). We use 5 x 10 s individual SPH particles to 
follow the gas dynamics. We model the supersonic turbu- 
lent motions that are observed to be present in molecular 
clouds by including a divergence-free random Gaussian ve- 
locity field with a power spectrum P(k) oc fc -4 where k 
is the wavenumber of the velocity perturbations (Ostriker, 
Stone & Gammie 2001). In three dimensions, this matches 
the observed variation with size of the velocity dispersion 
found in molecular clouds (Larson 1981). The velocities are 
normalised to make the kinetic energy equal to the abso- 
lute magnitude of the potential energy so that the cloud is 
marginally unbound. In contrast, the thermal energy is ini- 
tially only 1 per cent of the kinetic energy. Thus the cloud 
contains 1000 thermal Jeans masses, 

{5R 9 T\ 3/2 /4 \-^ 2 

where T is the gas temperature, R g is the gas constant, G is 
the gravitational constant, /i is the mean molecular weight 
and p is the density of the gas. Thus, in the absence of the 
turbulence and hence kinetic support, the cloud could be 
expected to fragment into of order 1000 stars should suffi- 
cient structure be present (Bonnell 1999). We force the gas 
to remain isothermal throughout the simulation, emulating 
the effect of efficient radiative cooling at the low gas den- 
sities present. We do not include any feedback (radiative 
or kinematic) from the newly formed stars. We expect that 
feedback, especially from massive stars, will start to become 



important by the end of the simulation. The simulation was 
carried out on the United Kingdom's Astrophysical Fluids 
Facility (UKAFF), a 128 CPU SGI Origin 3800 supercom- 
puter. 

Once fragmentation has produce a protostar, we replace 
the constituent SPH particles with a single sink-particle 
(Bate, Bonnell & Price 1995). These sink-particles, used to 
follow the newly formed stars, interact only through gravi- 
tational forces and by accretion of gas particles that fall into 
their sink-radii. Sink-particle creation occurs when the dens- 
est gas particle (at a given time) and its pb 50 neighbours 
are self-gravitating, subvirial and occupy a region smaller 
than the sink- radius (200 AU). For the simulation reported 
here, this requires a gas density of <T.5 x 10~ 15 g cm -3 . 
Gas particles that fall within a sink-radius of 200 AU are 
accreted if they are bound to the sink-particle whereas all 
gas particles that fall within 40 AU are accreted, regard- 
less of their properties. A minimum number of SPH par- 
ticles is required to resolve a fragmentation event (Bate & 
Burkert 1997) implying that our completeness limit for frag- 
mentation is approximately Q.IMq. Thus, we cannot resolve 
any protostars that would form with (initial) masses below 
this limit. Therefore, although we do resolve the bulk of the 
fragmentation, there will be a number of lower-mass stars 
and brown dwarfs that we do not capture in our simulation. 
Gas accretion onto the stars then increases their masses and 
removes SPH particles from the simulation. In order to min- 
imise computational expense, we smooth the gravitational 
forces between stars at distances of 160 AU. This means that 
only the widest binaries and multiple systems are resolved 
and that stellar collisions are not included in the simulation. 
The use of gravitational softenning allows us to evolve the 
system further than has previously been achieved and thus 
evaluate the formation process of the stellar cluster. 



3 HIERARCHICAL CLUSTER FORMATION 

The initial evolution of the molecular cloud is due to the 
turbulent motions present in the gas. The supersonic turbu- 
lence leads to the development of shocks in the gas, produc- 
ing filamentary structures (Bate, Bonnell & Bromm 2003). 
The shocks also remove kinetic energy (assumed to be ra- 
diated away) and thus remove the turbulent support locally 
(Ostriker et al. 2001). The chaotic nature of the turbulence 
leads to higher-density regions in the filamentary structures 
which, if they become self-gravitating, collapse to form stars. 

Star formation occurs simultaneously at several differ- 
ent locations in the cloud. Figure plots the column density 
through the molecular cloud at four different times over the 
4.5 x 10 5 years (2.6 free-fall times, tg) that we follow the 
evolution. The stars that form first are in the highest den- 
sity gas where the dynamical timescale is the shortest. This 
generally occurs in the deepest parts of local potential min- 
ima. Surrounding clumps with slightly lower gas densities 
form additional stars. Both the stars and the residual gas 
are attracted by their mutual gravitational forces and fall 
toward each other. Gas dynamics dampen the infall veloc- 
ities allowing the systems to rapidly merge to form a high 
density subcluster containing from five to several tens of 
stars. The number of stars in each subcluster increases as 
further star formation occurs nearby, and these stars fall 
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Figure 1. The stellar cluster forms through the hierarchical fragmentation of a turbulent molecular cloud. Each panel shows a region 1 
parsec on a side. The logarithm of the column density is plotted from a minimum of 0.025 (black) to a maximum of 250 (white) g cm -3 . 
The stars are indicated by the white dots. The four panels capture the evolution of the 1000 Mq system at times of 1.0, 1.4, 1.8 and 2.4 
initial free-fall times, where the free-fall time for the cloud is tg = 1.9 X 10 years. The turbulence causes shocks to form in the molecular 
cloud, dissipating kinetic energy and producing filamentary structure which fragment to form dense cores and individual stars (panel A). 
The stars fall towards local potential minima and hence form subclusters (panel B). These subclusters evolve by accreting more stars 
and gas, ejecting stars, and by mergers with other subclusters (panel C). The final state of the simulation is a single, centrally condensed 
cluster with little substructure (panel D). The cluster contains more than 400 stars and has a gas fraction of approximately 16 per cent. 



into the existing potential wells. This process repeats itself 
until several hundred stars are formed and mostly contained 
in five subclusters. The further evolution is marked by a de- 
creasing star formation rate as the subclusters, aided by the 
dissipative effects of their embedded gas, sink towards each 
other and finally merge to form one single cluster containing 
over 400 stars. The final cluster is approximately spherical 
in shape with a centrally condensed core as is observed in 
young stellar clusters (Hillenbrand 1997; Lada 1999). 



The hierarchical nature of the formation process is il- 
lustrated in figure |21 which shows the evolution of the local 
and global stellar number-densities for the cluster. The local 
stellar density is calculated for each star from the minimum 
volume required to contain the star's ten nearest neighbours. 
We use the median value of this distribution to quantify a 
typical local stellar density. In contrast, the global stellar 
number-density is calculated from the volume required to 
contain half of the total number of stars. This typifies the 
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Figure 2. The evolution of the local and global stellar number- 
densities are plotted as a function of time in years (the initial 
free-fall time is % = 1.9 X 10 5 years). The local stellar number- 
density (solid line) is calculated from the volume required to hold 
the stars' 10 nearest neighbours; the star having the median den- 
sity is shown. The global stellar number-density (dashed line) is 
calculated from the volume required to contain fully half of the 
total number of stars. This is equivalent to the local density for 
a monolithic formation process (see text). The rapid rise of the 
local (hierarchical) stellar density compared to the global (mono- 
lithic) stellar density is due to the formation of subclusters. The 
discrepancy between the two values decreases with time as the 
substructure is erased to produce a single cluster. 



stellar densities expected from a monolithic (homogeneous 
and structureless) formation scenario. We see that the local 
number-density increases rapidly once the first stars form 
and fall towards each other in their local subcluster. The 
local number density attains a maximum of 10 s stars pc -3 , 
up to 100 times that for a monolithic collapse. The differ- 
ence between the two values indicates that the stars occupy 
but a small fraction of the total volume of the star forming 
region. This has significant implications for the probability 
for interactions (see below). 

The local number-density decreases after reaching a 
maximum during the subclustering phase. This decrease is 
due to the ejection of stars from the subclusters through dy- 
namical interactions, and due to the kinetic heating during 
the merging of subclusters. This process erases the substruc- 
ture fairly quickly, producing a single, centrally condensed 
cluster. 



4 ACCRETION AND THE INITIAL MASS 
FUNCTION 

In addition to acting as a reservoir for star formation 
and as a damping force of the stellar dynamics, gas is ac- 
creted onto individual stars, thereby increasing their masses. 
The stars compete for the gas, with those in the bottom of 



Figure 3. The final mass distribution of stars is plotted as a 
function of the logarithm of the mass. The distribution uses bins 
in the logarithm of the mass such that a Salpeter IMF has a slope 
of r = —1.35. The minimum mass for the simulation is plotted 
as the dashed line while the diagonal line notes a T = — 1 slope. 
The higher-mass stars appear to have a steeper distribution while 
intermediate-mass stars have a shallower one. Stars below w lMg 
have an approximately flat distribution. 



their local potential wells accreting the most and becoming 
the most massive stars (Bonnell et al. 1997, 2001a). Thus, 
the first stars to form are frequently the most massive due 
to this accretion process. This process, termed competitive 
accretion, is a leading candidate to explain the apparently 
universal initial mass function (IMF) of stars (Bonnell et 
al. 2001b; Klessen 2001). 

Here, gas accretion results in final stellar masses that 
range from approximately 0.07 to 27 Mq . The final mass dis- 
tribution, shown in ngure[3]is consistent with observed IMFs 
(Hillenbrand 1997; Luhman et al. 2000; Meyer et al. 2000). 
The distribution has a near flat slope for low-mass stars, 
that turns into an increasingly steeper slope for more mas- 
sive stars (see also Bonnell et al. 2001b; Bate et al. 2003). 
The higher-mass distribution is broadly consistent with a 
r ~ — 1 slope (dN(logm) = m r dm, where the Salpeter 
slope has F = —1.35) although could also be fit by a shal- 
lower slope for intermediate-mass stars and a steeper slope 
for high-mass stars. This high-mass slope is similar to the 
recent result where high-mass stars are formed through a 
combination of gas accretion and stellar mergers (Bonnell 
& Bate 2002). In this simulation, the gravitational poten- 
tial wass not softened and binary systems, formed through 
three-body capture, had separations as small as 10 AU. 

The median and mean stellar mass are 0.43 and 1.38 
Mq, respectively. At the end of the simulation, 42 per cent 
of the total mass remains in gas, although much of the gas 
is no longer bound to the cluster due to the initial turbu- 
lence. The cluster contains 494 Mq in a 0.25 pc radius, but 
only 16 per cent is in the form of gas. This star formation 
efficiency, although comparable to that observed for young 



Hierarchical cluster formation 5 



stellar clusters (Lada 1991; Lada & Lada 1995; Clarke et 
al. 2000), does not include any background molecular cloud 
not directly involved in the cluster formation process. Fur- 
thermore, as the fraction of mass in stars is an ever increas- 
ing function of time, its value depends on when we halt the 
simulation. At some point in the process, feedback from the 
more massive stars is expected to clear the remaining gas 
from the system. This gas expulsion will force the cluster to 
expand as may be occurring in the ONC (Kroupa, Hurley & 
Aarseth 2001). Feedback from massive stars could also af- 
fect the accretion process. As feedback acts relatively quickly 
and only once the massive stars have formed, its most prob- 
able effect will be a freezing of the resultant mass function 
(Bonnell et al. 2001b). 



5 DISCUSSION 

The hierarchical nature of the formation process has many 
interesting implications for star formation. Subclustering 
means that individual stars are in regions of higher stellar 
number-density than they would be for a monolithic forma- 
tion process (Fig.H Fall & Rees 1985; Scally & Clarke 2002). 
The high number- density of stars, coupled with the rela- 
tively small number of stars in each subcluster, and thus 
smaller velocity dispersion, results in closer, and stronger 
stellar interactions than would otherwise occur (Scally & 
Clarke 2002). Such stellar interactions can harden binaries to 
explain the closest systems (Bate, Bonnell & Bromm 2002), 
truncate circumstellar discs (Hall, Clarke & Pringle 1996; 
Bate et al. 2003) decreasing their masses and thus lifetimes, 
trigger fragmentation in the disc (Boffin et al. 1998; Watkins 
et al. 1998) and possibly even liberate a population of plan- 
ets from their parent stars if planets form quickly (Bonnell 
et al. 2001c; Hurley & Shara 2002). The maximum number- 
density of stars is sufficiently high (10 7 to 10 s stars pc -3 ) 
that stellar mergers may play a role in forming the most 
massive stars (Bonnell, Bate & Zinnecker 1998; Bonnell & 
Bate 2002). 

Figure [I] plots the distribution of closest approaches for 
each of the 418 stars formed in the simulation. This distribu- 
tion is calculated, for each star, as the minimum distance to 
any other passing star sometime during the evolution. The 
distribution extends from 10 AU to > 10 4 AU. The small 
peak at large separations indicates the few stars that form 
in relative isolation in the molecular cloud, and never enter 
into a cluster. Nearly half of the stars in the simulation have 
interactions within the 160 AU resolution limit where we 
start to smooth the gravitational forces. These are therefore 
upper limits for the actual closest approaches. 

We see in figure 01 that approximately one third 
(140/400) of the stars have encountered another star within 
100 AU, sufficiently close to perturb the circumstellar disc 
(McCaughrean & O'Dell 1996) or a binary system (Duquen- 
noy & Mayor 1991). A close passage of a star is likely to 
disrupt a circumstellar disc down to one third of the mini- 
mum separation between stars (Hall et al. 1996) , thus limit- 
ing their mass reservoirs and lifetimes (although subsequent 
accretion may replenish the discs [Bate et al. 2003]). Fig- 
ure 0] also plots the corresponding distribution of closest ap- 
proaches for more massive stars, with m > 3Mq. The great 
majority (34/40) of these massive stars have had a close in- 
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Figure 4. The distribution of minimum closest approaches for 
all stars (solid line) and stars more than three Mq (dashed line) is 
plotted as a function of the logarithm of the minimum separation. 
Gravitational softenning occurs for approaches within 160 AU. 
Thus, the measured values for the closest approaches are upper 
limits in these cases. 



teraction within 100 AU. Interactions within 100 AU result 
in a disc truncated to 30 AU or less (Hall et al. 1996). This 
suggests that many of the stars, and virtually all the mas- 
sive stars, found in young stellar clusters, should have small 
(less than 30 AU) or non-existent discs. A similar result was 
found in the smaller cluster simulation of Bate et al. (2003) 
which resolved discs down to ~ 10 AU. 



6 CONCLUSIONS 

The results of the numerical simulation presented here point 
to a hierarchical formation process for stellar clusters. Star 
formation occurs along filamentary structures engendered 
in the molecular cloud due to supersonic turbulent mo- 
tions. The stars fall together to form many small subclus- 
ters. These subclusters grow by accreting other stars and 
eventually merge to form a large-scale cluster containing 
over 400 stars. The heirarchical nature of the formation pro- 
cess means that close interactions between forming stars is 
much more important than would be in a monolithic for- 
mation process. Thus, stars are more likely to have suffered 
close interactions that can truncate their circumstellar discs, 
harden existing binaries and possibly liberate planets from 
their parent stars. In particular, the hierarchical process and 
the corresponding higher stellar number-densities imply that 
approximately one third of low-mass stars, and most high 
mass stars, should have had their discs truncated to within 
30 AU. 

Evidence for a hierarchical formation process exists in 
the results of Testi et al.(2000) showing the hierarchical or- 
ganisation of the sub-mm cores in the Serpens molecular 
cloud. Detection of the subclustered phase at optical and 
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near infrared wavelengths is problematic due to the large 
dust extinction present in the cloud. Fortunately, the up- 
coming Space Infrared Telescope Facility (SIRTF) will ob- 
serve star forming regions in the mid-infrared and should 
be able to determine the level of hierarchical subclustering 
present in even the youngest stellar clusters. 

Finally, this numerical simulation demonstrates that 
there is a strong link between the large-scale dynamics of 
star formation and the small-scale stellar and protoplane- 
tary disc properties. This highlights the importance of a 
global approach to the problem. Future studies will need to 
include even larger-scale processes such as the formation and 
evolution of molecular clouds. In this way, we will be able 
to determine how and when star formation is initiated and 
thus tie star formation into Galactic evolution. 
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